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Abstract This paper provides an improved flamelet/progress variable (FPV) model for 
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Z, and of the progress parameter, A. Steady-state FPV models are built presuming the func¬ 
tional shape of the joint PDF of Z and A in order to evaluate Favre-averages of thermody¬ 
namic quantities. The mixture fraction is widely assumed to behave as a passive scalar with 
a mono-modal behaviour modelled by a j8-distribution. Moreover, under the hypothesis that 
Z and A are statistically independent, the joint PDF coincides with the product of the two 
marginal PDFs. In this work we discuss these two constitutive hypotheses. The proposed 
model evaluates the most probable joint distribution of Z and A, relaxing some crucial as¬ 
sumption on their statistical behaviour. This provides a more general model in the context of 
FPV approach and an effective tool to verify the adequateness of widely used hypotheses. 
The model is validated versus experimental data of well-known test cases, namely, the San- 
dia flames. The results are also compared with those obtained by the standard FPV approach, 
analysing the role of the PDF functional form on turbulent combustion simulations. 

Keywords Presumed PDF • Methane combustion • Turbulent non-premixed combustion • 
Reynolds-Averaged Navier-Stokes equations 


1 Introduction 

The development of new technologies to enhance and control combustion processes is nowa¬ 
days fundamental in order to increase efficiency and reduce emissions in many engineering 
applications, such as internal combustion engines, advanced gas turbine systems, pulse det¬ 
onation engines, high-speed air-breathing propulsion devices. In this context, the simulation 
of turbulent reacting flows is very useful to reduce experimental costs and to advance the 
comprehension of the basic physical mechanisms. Turbulent combustion is a multi-scale 
problem, where the interaction between chemical kinetics, molecular, and turbulent trans¬ 
port occurs over a wide range of length and time scales. The numerical simulation of such 
phenomena with detailed chemistry is today still prohibitive. The huge computational cost 
stems from the fact that, even for a simple fuels, detailed kinetic mechanism can involve 
thousands of reactions and, consequently, hundreds of chenncal species. One way to over¬ 
come this problem is to simplify it using reduced chemical models with a limited number of 
reactions and species involved, but still the number of degrees of freedom may be prohibitive 
for applications of engineering interest. Therefore, several approaches have been proposed 
that simplify the problem by pre-computing, in terms of a small number of variables, the 
entire chemical state-space mi2iia4i . Among such models, the fiamelet-progress variable 
(FPV) approach tan is considered in the present work. For the case of non-premixed com¬ 
bustion of interest here, the basic assumptions of the fiamelet model are fulfilled for suffi¬ 
ciently large Damkdhler number. Da. In fact, when the reaction zone thickness is very thin 
with respect to the Kolmogorov length scale, turbulent structures are unable to penetrate into 
the reaction zone and cannot destroy the laminar flame structure. Effects of turbulence only 
result in a deformation and straining of the flame sheet and locally the flame structure can 
be described as function of the mixture fraction, Z, the scalar dissipation rate, the 

time. The scalar dissipation rate, % — 2Dz(VZ)^, is a measure of the gradient of the mixture 
fraction representing the molecular diffusion of the species in the flame region, Dz being the 
molecular diffusion coefficient. Therefore, the entire flame behaviour can be obtained as a 
combination of solutions of the laminar fiamelet equations. In the present work we consider 
a further simplification assuming a steady fiamelet behaviour, so that chemical effects are 
entirely determined by the value of Z. Moreover, % defines the effects of the flow on the 
flame structure according to the following steady laminar fiamelet equation (SLFE) for 
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the generic variable 0 


^ 2 dz^ 


= , 


( 1 ) 


where p is the density and 6)^ is the source term related to 0 IJ), different from zero in the 
case of finite rate chemistry. In this work, we employ the FPV model proposed by Pierce 
and Moin f/llSl to evaluate all of the thermo-chemical quantities involved in the combus¬ 
tion process. This approach is based on the parametrization of the generic thermo-chemical 
quantity, (j), in terms of the mixture fraction, Z, and of the progress parameter. A, instead of 


Z: 


0=F>(Z,A). 


( 2 ) 


Using such a parameter, independent of the mixture fraction, one can uniquely identify 
each flame state along the stable and unstable branches of the S-shaped curve. A suitable 
definition of A leads to a substantial simplification of the presumed PDF-closure model. In 
contrast, the solution of the transport equation for A is quite complex since it requires non¬ 
trivial modelling of several unclosed terms |8( . In order to overcome such a difficulty, the 
progress parameter is derived from a reaction progress variable, C, such as the temperature 
or a linear combination of the main reaction products, whose behaviour is governed by a 
simpler transport equation. Therefore, a transport equation for C is solved and the flamelet 
library is parametrized in terms of Z and C. Requiring that the transformation between A 
and C be bijective, from equation O one has 


A=FcHz,C), 


(3) 


and any thermo-chemical variable can then be expressed as: 




(4) 


The choice of the progress variable is not unique and some recent works discuss in details 
this issue proposing a procedure for the selection of the optimal definition of C iQllTolfTTl . 
Here the progress variable is defined by the sum of the mass fractions of the main products; 
for methane combustion: 

C = Yh,o + Yco + Yco2 • ( 5 ) 


This turns out to be a suitable choice for the Sandia flames test case considered in this 
work 121 . Observe that, since only the main products of a combustion process are considered 
the maximum value assumed by the progress variable is smaller than unity. A stretching with 
respect to the minimum and maximum conditioned value of C over the mixture fraction is 
made, namely: 


A = 


C CMin\'Y^ 
^Max\Y‘ CmIu \Yi 


( 6 ) 


This definition greatly simplifies the presumed PDF closure procedure, recovering the uni¬ 
tary length of the variation interval for the progress parameter fT^ and so overcoming 
the issue in the integrals evaluation. Figure shows the effect of such a normalization. 
Equation 0 is taken as the solution of the SLFE 0. Each solution corresponding to a 
given value of z is a flamelet and the solution variety over x — Xst is shown in figure 
From equation ^ one can obtain the Favre-averaged value of (j) and of its variance using 
the definitions: 

^=J J F^{Z,A)P{Z,A)dZdA, 


( 7 ) 
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Fig. 1 Schematic of the procedure for the conditional normalization of the progress variable. 



Fig. 2 Distribution obtained plotting the maximum temperature of each flamelet versus the scalar stoichio¬ 
metric dissipation rate of about 250 flamelet with locally reflned distribution in Xst direction. 


= J J[F^{Z,A)-^f P(Z,A)dZdA, 


( 8 ) 


where P(Z, A) is the density-weighted PDF. This function plays a crucial role in the defini¬ 
tion of the model, affecting both its accuracy and computational cost. Moreover, the choice 
of P(Z,A) is not straightforward because of the unknown statistical behaviour of the two 
variables Z and A and its definition is still an open problem whose solution is being 
pursued by several researches ISHSlllSIlldlfTSl . 

The aim of this work is to provide an analysis of state-of-the-art FPV models and to pro¬ 
pose a more general scheme based on the statistically most likely distribution (SMLD) ca 
approach for the joint PDF of Z and A. In fact, this allows one to avoid any assumption about 
the statistical correlation of the two variables and about the behaviour of the mixture frac¬ 
tion. In order to assess the capability of the combustion models, they are employed here in 
conjunction with a Reynolds-Averaged Navier-Stokes (RANS) solver with a k-co turbulence 
model closure. 

This work is organised as follow: in section [T] a brief introduction of the state-of-the- 
art is presented in order to show the motivation and the basic outlines of flamelet models; 
then, in sections and the theoretical and numerical details of the models employed are 
provided. The comparison between four models for the definition of the presumed-PDF is 
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discussed in sectionj^and numerical results obtained studying two subsonic flames, namely 
Sandia Laboratories flame D and E, are presented. Finally, some conclusions are provided. 


2 The flamelet progress variable models 

2.1 The standard FPV model (model A) 


In this section, the standard FPV model (called here model A) is briefly described. It is based 
on the following hypotheses (Hp). 

Hp \ \ The mixture fraction Z and the progress parameter A are the two independent vari¬ 
ables of the combustion model; they form a basis by which one can derive all of the thermo¬ 
chemical quantities. 

Hp 2\ The steady flamelet assumption is used, so that the solution variety of equation 0 
is made over % — Xst, see figure A key difficulty to integrate the flamelet equations is 
to know a priori information about the scalar dissipation rate dependence on the mixture 
fraction x{Z). In fact, flamelet libraries are computed in advance and are assumed to be 
independent of the flow field. Therefore, the dependence of on Z has to be modelled El 
[TslfT^ . In this work, the functional form of ;f(Z) has been taken from an idealized con- 
figuration. Peters EOl developed an analytic expression for the distribution of the scalar 
dissipation rate in a counterflow diffusion flame, which leads to an inverse error functional; 


X{Z)=Xst 


<P{Zs,) ’ 


(9) 


where Zgt is the mixture fraction evaluated at the stoichiometric point and ^(Z) is the dis¬ 
tribution of scalar dissipation rate for the counterflow diffusion flame |[20l. 

Hp 3 : The PDF in Q and ^ is presumed and its choice establishes the statistical correlation 
between Z and A. Therefore, employing Bayes’ theorem. 


P(Z,A)=P(Z)P(A|Z), 


( 10 ) 


to evaluate such a PDF, one usually presumes the functional shape of the marginal PDF of 
Z and of the conditional PDF of A |Z. 

Hp 4: The statistical behaviour of the mixture fraction is described by a j8-distribution. In 
fact, even though the definition of P(Z) is still an open question (Tb), it has been shown by 
several authors that the mixture fraction behaves like a conserved scalar whose statistical 
distribution can be approximated by a j8-function 1121112211231 . The two parameter family of 
the j8-distribution in the interval x G [0,1] is given by: 


j8(x;V,V'2)=x^-i(l-x)^-i 


r{a + h) 

r{a)r{hy 


( 11 ) 


where F (x) is the Gamma-function and a and b are two parameters related to x and x"^ 


v //2 


v //2 


Hp 5: Statistical independence of Z and A is assumed, so that P{Z,A) — P{Z)P{A). 

Hp 6: P(A) is a Dirac distribution, implying a great simplification in the theoretical frame¬ 
work. 
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Under these assumptions, it can be shown that there is only one solution of equation 
0 for each value of the scalar dissipation rate. With these criteria, the Favre-average of a 
generic thermo-chemical quantity is given by: 


(l> = j j F^{Z,A)^{Z)8{A-A)dZdA = J F^{Z,A)l5iZ)dZ. (13) 


Therefore, one has only three additional transport equations (for Z, and C) to evaluate 
all of the thermo-chemical quantities in the flow, thus avoiding the expensive solution of one 
transport equation for each chemical species. 

2.2 The FPV model with j8-distribution for P(A) (model B) 

It is well known that a reactive scalar, such as C, depends on a combination of solutions 
of equation 0 for each chemical state and therefore its PDF cannot be accurately ap¬ 
proximated by a Dirac distribution lH). Thus, a different model (called here model B) has 
been proposed in l24l . in which Z and A behave in the same way, namely, following a 
j8-distribution: 


P(Z,A) = /3(Z)/3(A). 


(14) 


This reformulation of Hp 6 can be considered as an improvement with respect to model A. 
However, this assumption has an additional cost since model B requires the evaluation of an 
additional transport equation for C"'^. 

2.3 The FPV model with SMLD for P(A) (model C) 

The probability distribution of a reacting scalar is often multi-modal and its functional form 
depends on the turbulence-chemistry interaction. It is noteworthy that the j8-function is suit¬ 
able to describe the behaviour of a passive scalar, such as Z; in fact, it becomes bi-modal for 
sufficiently large values of the variance. In order to have an improved evaluation of P(A), 
a more general distribution that allows multi-modal behaviour even for a small variance is 
needed. Therefore, Ihme and Pitsch P8( proposed a new model evaluating the conditional 
PDF as the statistically most likely distribution (SMLD) based on the values of Z, Z"^, A 
and A"^. The rational behind the SMLD approach is that one can construct a PDF based on 
the knowledge of a small number of moments and on the minimization of the uncertainty 
in order to have the most probable distribution. In this perspective, model C is more general 
than model A and model B, since the functional shape of P(A) is not fixed. Relying on the 
knowledge of the first three moments of P(A), a unique measure, S, of the predictability of 
a thermodynamic state can be defined II16II25I . S is an entropy function depending on P(A), 
S — S{P{A)) f26i , that can be thought of as the Boltzmann’s entropy: 



(15) 


where 2(A) is a density function proposed by Pope lITbl that takes into account the bias in 
the composition space when no information about the moments are given. 



( 16 ) 
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The goal is to construct a PDF that maximizes the entropy S. Following the Lagrangian opti¬ 
mization approach, the functional 5* is defined by involving the constraints on the moments: 







(17) 


In the above equation fin are the Lagrange’s multipliers while the last fraction term is intro¬ 
duced to normalize P(A). Equation ( p^ can be obtained from equation GD independently 
of the particular form of 2(A) considering the analytical form of the mean and the variance 
as constraints. The term depending on fin corresponds to the two constraints given by the 
first and the second moments of the PDF, namely A and A"^. The first moment is enforced 
considering the last term, needed to ensure the unity sum of P(A). The expression for 
P(A), obtained by maximising 5*, reads: 

P(A) = e(A)Lexp|-£^(A-;vr|, ( 18 ) 

Mo ^ J 

where: 


Mo = 

1 dAP{A), 

Jo 

(19) 

-Ml = 

['dAdA{PiA))=P{l)-P{0), 

Jo 

(20) 

M2 A = 

f dAdA[{A-A)P{A)]=P{l)-Anu 

Jo 

(21) 


since A is bounded in [0,1]. The model still needs an additional assumption for closure. 
Here, we make a different choice with respect to the original SMLD model of Ihme and 
Pitsch m. We assume that the first and the last point of P(A) are equal to the first and last 
points of j8 (A) evaluated with the same values of the mean and variance: 


P(1;A,A"2) = j8(l;A,A"2), P(0;A,A"2) = j8(0;A,A"2). (22) 


This assumption does not affect the multi-modal nature of the distribution, but simpli¬ 
fies the model implementation since there is no need to evaluate the roots of a non-linear 
system at each computational point as in the original version of the SMLD model llsl . 
Figu re[3] shows the comparison between Psml,2, evaluated with the constraint given by equa¬ 
tion ( |2^ (model C), and the distribution obtained by solving the natural constraint given by 
the non-linear system dij — fi2ik^kj = 0 HI, where is the covariance matrix, namely 
PsML 2- distributions are practically coincident for low and high values of the vari¬ 

ance. On the other hand, for its intermediate values, using the constraints in equation ( [^ 
provides a different behaviour of the two PDFs at the extrema of the range of variation of Z 
and A. However, such a discrepancy is localized at the extrema and only slightly affects the 
distribution for intermediate values, so that the simplified constraints in equation ( [22] ) can 
be considered satisfactory. 







A. Coclite et al. 


(a) Z=0.3 Z"^=0.1 



(c) Z=0.3 Z"^=0.01 



(b) A=0.8 A"^=0.04 




Fig. 3 Comparison among three PDFs for two set of values of the mean and the variance: j8 distribution 
(red dashed line); Psml,2 (black solid line) and P^mli El (symbols). Boundary values: Panel a, Psml,2{^) = 
i3(0) =40 ,Psml. 2(1) =i3(l) = 17,Pi^i_2(0) = =0; Panel b.PsML, 2 ( 0 ) = i3(0) = 0,P5Mi.2(l) = 

i3(l) = 13, = 0,P^i_ 2 (l) = 0; Panelc, PsMik^) =m = 0. PsMLaW = i3(l) = 0,P'Si_2(0) = 0, 

^™1,2(1) = 0; Panel d, PsmlM = m= 0, Psml.iW =m = 0, = 0, PiJ?i, 2 (l) = 1-13. 


2.4 The FPV model with SMLD for P(Z, A) (model D) 


The major advantage of the SMLD approach over usually employed presumed PDF closure 
models is that it provides a systematic framework to incorporate an arbitrary number of mo¬ 
ment information. It has been shown jS) that, even though equation 0 is based on assuming 
as independent variables Z and A, one has to properly take into account the statistical cor¬ 
relation between Z and A due to the effect of turbulence modelling, employing a suitable 
estimation of the joint PDF. Therefore, we propose a different model (model D) in which the 
SMLD approach is applied directly to the joint distribution. In this way one does not need 
hypotheses 4-6 and the most probable PDF can be evaluated considering only the known in¬ 
formations about the statistical means of the joint distribution. Assuming that the first three 
moments of the joint probability density function P(x), where x = (Z, A)^, are known, the 
procedure employed for model C is here extended to the case of the two dimensional PDF, 
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obtaining: 


PsML ,2 (Z, A) = exp{ - [mi,i (Z - Z) + ^ 1 , 2 (A - A)] 

- \ [m2,ii(Z-Z)2 + M2,12(Z-Z)(A -A) 

+ M 2,21 (A - A) (Z - Z) + M 2,22 (A - A )2] } . (23) 

In the equation above, Q{A\Z) — Q{A) is given by equation being the bias PDF for 
Z a constant (assumed equal to one), as discussed by Pope fTJ] due to the linearity of the 
mixture fraction transport equation in Z. Moreover , jUq is a scalar, /ii is a two - component 
vector, and ^ is a square matrix of rank two: 

Mo = y dxPsMiai'!^), (24) 

“Ml,; = j dxdxiPsML,2{^) 

= (25) 

4/ - M 2 ,in = Jdxdx^{{xi - ^l)PsML, 2 {x)) 

= /3(l;|,,|^')-|iMi./, (26) 

where /, k, n, and / indicate the vector components; (^/ and are the mean, the fluctua¬ 
tion and the variance of the i-ih component of x, respectively v/ — and 
To proof equations <24}, ( [^ and ( [^ one needs to consider Bayes’ theorem applied to 

PsML.li^) = PsML,2{m)PsML,2{^2\m) , (27) 

where P 5 ml, 2 (-^i) is the marginal PDF and P 5 ml ,2 (-^2 l-^i) is the conditional one. These two 
functions are both PDFs, therefore by definition one has: 


J dxiPsML,2{2Cl) — 1, (28) 

J dX2PsML,2{2C2\2Cl) = T (29) 

Considering a probability space (^2 ,F,P), where Q is the sample space, F the events algebra 
and P the PDF, and two events A,B e Q with P(A) ^ 0) the conditional probability of 
5 at a given A is: P{B\A). P(-|A) : F ^ [0,1] is actually a PDF as proved by invoking the 
characterization theorem, so that, P{Q |A) = 1. Therefore, (^2, F, P(- |A)) is a new probability 
space in which there is the information that A is an event always verified. In some way one 
can think that A is a sort of parameter that measures the known region of Q . 

Let us consider equation ( [25] ): 

-Mi,i = J dxdxiPsML,2ix) = 

= PsML, 2 {xi) lo • 


(30) 
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Fig. 4 Comparison between j8(Z) and PsuL.ii^), and j8(A|Z) and PsML.ii^l^) distributions at two fixed 
arbitrary mean values, namely Z = 0.5 and A |Z = 0.8 value and for three values of the two variances, respec¬ 
tively. 


Straightforwardly comes the proof of the equations ( [^ and trivially that of the equa¬ 
tion 

It is noteworthy that, assuming the statistical independence of Z and A and a j8-distribution 
for P(Z), model D automatically reduces to model C. Indeed, in such a case. 


PsmmiZ.A) = P{Z)PsmlM^) = P{Z)Psml,2(^) , (31) 

where the first multiplier jUq is still given by equation ( [^ , while the second and the third 
ones, /ii and reduce to scalar quantities given by equations ( [20] ) and dB- respectively. 
Figure |4| shows the comparison between the two PDFs for (arbitrarily chosen) Z = 0.5 (left 
panel) and A |Z = 0.8 (right panel). Three curves are plotted in each panel, corresponding to 
three values of the variance in the interval [0,Z(1 — Z)] = [0,0.25] and [0,A|Z(1 — A|Z)] = 
[0,0.16], respectively. One can see that the employed closure for the Lagrangian multiplier 
at the extrema (see equation ([^), does not prevent the SML distribution to separate signifi¬ 
cantly from the j8-distribution for high values of the variance. In fact, the two PDFs slightly 
differ for small and large value of the variance, the major differences being obtained for 
value not so close to the upper and lower limits of the admissible interval. 


3 Governing equations 

3.1 Flow equations and numerical solution procedure 

The numerical method developed in ll28l has been employed to solve the steady-state RANS 
equations with the k-(0 turbulence closure. For multi-component reacting compressible flows, 
the system of governing equations can be written as: 


a,Q + V-(E-Ev) = S, (32) 

where t is the time. This equation is solved in an axisymmetric cylindrical reference frame, x 
and y representing the axial and the radial coordinate, respectively; Q=(p,pUx,puy,pH,pk, 
pco,pRn) is the vector of the conserved variables; p, {ux, Uy), H indicate the Favre-averaged 
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values of density, velocity components and specific total enthalpy, respectively; k and (O are 
the turbulence kinetic energy and its specific dissipation rate; Rn is a generic set of conserved 
variables related to the combustion model. In this framework, is the set of independent 
variables of the flamelet model, namely, Z, Z"^, C, E and Ey are the inviscid and vis¬ 
cous flux vectors 1291 , respectively; S is the vector of the source terms. 

A cell-centred finite volume space discretization is used on a multi-block structured mesh. 
The convective and viscous terms are discretized by the third-order-accurate Steger and 
Warming 1^ flux-vector-splitting scheme and by second-order-accurate central differences, 
respectively. An implicit time marching procedure is used with a factorization based on the 
diagonalization procedure of Pulliamm and Chaus see ED, SO as to allow a standard scalar 
alternating direction implicit (ADI) solution procedure t37\ . Only steady flows are dealt with 
in this paper. Therefore, the unsteady terms are dropped from the governing equations and 
the ADI scheme is iterated in the pseudo-time until a residual drop of at least five orders of 
magnitude for all of the conservation-law equations ( [^ is achieved. Characteristic bound¬ 
ary conditions for the flow variables are imposed at inflow and outflow points, whereas no 
slip and adiabatic conditions are imposed at walls; k, ft), and Rn are assigned at inflow points, 
whilst they are linearly extrapolated at outflow points. At solid walls, k is set to zero and ft) 
is evaluated as proposed by |33l: 


m = 60 


V 

0 . 093'2 , ’ 


(33) 


where is the distance of the first cell center from the wall; the homogeneous Neumann 
boundary condition is used for Rn (non-catalytic wall). Finally, symmetry conditions are 
imposed at the axis. 


3.2 Turbulent FPV transport equations 

For the case of a turbulent flame, the solution of the SLFE, namely equation D, must be 
written in terms of the Favre averages of Z and C and of their variance. The mean and the 
variance of the progress parameter are computed using equation ([^. Using model A one 
can tabulate all chemical quantities in terms of Z, Z"^ and C because of the properties of 
the 5-distribution. On the other hand, models B, C and D express (j) in terms of C"^ too and 
therefore they need to solve a transport equation also for C"^. The transport equations read: 


a,(pz) + v-(^z) = v- [(Z)+Z)ypvz], 

(34) 

5,(pZ"2) + V-(pSF2) = V [(£> + £>^)pVZ"2] - 


- pz + 2p4(VZ)2, 

(35) 

dt{pC) + V • (puC) = V • |^(Z) + Z)~)pVcj + pwc, 

(36) 

a,(pC^) + V.(puC^) = V. [(7) + Z)^)pVC^] - 


- p^ + 2p£>gVC)2 + 2pC^", 

(37) 


where Xc is modelled in terms of Z"^ and C"'^ [TS), namely Xc — ^ i^ diffusion 

coefficient for all of the species, given as D = v/Pr evaluated assuming a unity Lewis num¬ 
ber; V and Pr are the kinematic viscosity and the Prandtl number of the fluid, respectively; 
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Fig. 5 Detail of the computational grid close to the injector (upper frame) and overall view of the computa¬ 
tional domain (lower frame). 


Z)~ = D~ — D~ — D~ — Vf/Scf are the turbulent mass diffusion coefficients and Sct the 

Z ^"2 C C"2 ' 

Schmidt turbulent number setted equal to 0.9; (be is the source term for the progress vari¬ 
able precomputed and tabulated in the flamelet library. The gradient transport assumption 
for turbulent fluxes is used and the mean scalar dissipation rate, X and Xc, appear as a sink 
term in equations ( [^ and ( [?7] ), respectively. 

At every iteration, the values of the flamelet variables are updated using equations ( [34] )- 
and the Favre-averaged thermo-chemical quantities are computed, using equation (^. Such 
solutions provide the mean mass fractions which are used to evaluate all of the transport 
properties of the fluid and then to integrate the flow equations ( [^ . 


4 Numerical results 

This section provides the comparison among the results obtained using the four combus¬ 
tion models in order to assess the influence of the choice of the PDF on the prediction of 
turbulent non-premixed flames. The well-known Sandia flames are considered as suitable 
test cases, whose experimental data are available in the literature 1341 . The steady flamelet 
evaluations have been performed using the FlameMaster code Ea. 


4.1 Sandia Flames D and E 

In the Sandia flame experiments, piloted partially premixed CH 4 /air diffusion flames at the 
same pressure, equal to 100.6 kPa, and at different Reynolds numbers have been studied 
by Barlow & Frank 13611341 . The Reynolds number. Re, is based on the nozzle diameter, 
the jet bulk velocity, and the kinematic viscosity of the fuel. The diameter of the nozzle 
of the central jet is dref = 7.2 mm and the inner and outer diameters of the annular pi¬ 
lot nozzle are equal to 7.7 mm and 18.2 mm, respectively. The oxidizer air (7^2 = 0.233, 
¥^2 — 0.767) is supplied as a co-flow at 291 K. The main fuel jet consists of a mixture of 
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Fig. 6 Flame D thermo-chemical distributions along the axis {y/dref = 0). Model A, red dashed line; 
Model B, dashed-dotted line; Model C, dashed-dotted-dotted line; Model D, solid line; Symbols, experi¬ 
mental data. 



Fig. 7 Flame E thermo-chemical distributions along the axis {y/dref = 0)- Model A, red dashed line; 
Model B, dashed-dotted line; Model C, dashed-dotted-dotted line; Model D, solid line; Symbols, experi¬ 
mental data. 
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Table 1 Conditions for Sandia D-E flame experiment. 



u (m/s) 

7(K) 

Z 

Z"2 

Tu/„ 

(mm) 

Main jet 

49.6 ± 2 (D) 
74.4 ± 2 (E) 

294.0 

1 

0 

5% 

0.5 

Pilot stream 

11.4 ±0.5 (D) 
17.1 ±0.75 (E) 

1800 

0.27 

0.0075 

5% 

0.5 

Co-flow 

0.9 ± 0.05 

291.0 

0 

0 

0.5% 

1.0 


methane and air with a volumetric ratio equal to 1:3, and a stoichiometric mixture frac¬ 
tion Zst — 0.351; the corresponding mass fractions are Ych/^ — 0.156, 7^2 = 0.196, and 
Fa ^2 = 0.648. The pilot stream is a lean premixed gas mixture with an equivalence ratio 
of 0 = 0.77, corresponding to the equilibrium composition with Yco 2 =0.110, 7^2 = 0.056, 
7a^2 = 0.734, Yqh — 0.002, and YH 2 O — 0.098. The inlet value of the progress variable is 
equal to 0.208. Flame D {Re = 22400) presents very low degree of local extinction, whereas 
Flame E {Re — 33600) has significant and increasing probability of local extinction near 
the pilot. The experimental values of the axial velocities and temperatures have been as¬ 
signed at the inlet points (left boundary) together with the turbulent intensity levels, Tu/„, 
and turbulent characteristic length scale, as reported in Table Free-slip conditions 
have been imposed at the points on the upper boundary. The computational domain, shown 
in figure]^ is axisymmetric and includes a part of the burner; it has a length of 150 dref and 
27 dref along the axial and radial directions, respectively, and has been discretized using 
about 47000 cells. Computations have been performed using the combustion scheme de¬ 
scribed by the GRIMECH 3.0 t37l: 325 sub-reactions upon 53 specif. Thejlamelet library 
is computed over a grid with 125 uniformly distributed points in the Z and C directions and 
25 uniformly distributed points in the Z"^ and C'^ directions. 

Experimental data for the time-averaged major species and temperature are available 
as radial distributions at several stream-wise locations and as axial distributions along the 
center-line 03 The statistical errors of the measurements are below 10% for the CO mass 
fraction, and below 5% for the other quantities. Figures[^and[^show the comparison among 
the distributions of 7, Z, Z"^, and YH 2 O obtained with the four FPV models along the flame 
center-line for flame D and flame E, respectively. All models provide temperature distribu¬ 
tions in reasonable agreement with the experimental data, model A predicting a peak value 
of the temperature remarkably higher than the others. Moreover, for flame E a significant 
shift of the location of temperature peak is obtained with all models, probably due to the 
higher inlet velocity with respect to the flame D. The results of model B and model C are 
almost always overlapped. On the other hand, model D shows a slight improvement with 
respect to the others. This two figures show that all of the models are able to provide a 
good estimation of the mixture fraction distribution. Moreover the evaluation of the H 2 O 
mass fraction is in good agreement with the experimental data and here the four models 
present very small differences since they are based on the same efficient kinetic mecha¬ 
nism GD. The distributions of temperature and mixture fraction obtained with model C and 
model D at three axial sections are shown in figure for flame D and E. The comparison 
with model A is provided only for the first section, where a large discrepancy is observed in 
the temperature profile. In fact, one can see that model A predicts a maximum temperature 
of about 2200 K at x/dref — 1 whereas the maximum experimental value is about 1900 K\ 
in contrast, models C and D provide an improved predicted value of about 2000 K. Mod- 






An SMLD joint PDF model for turbulent non-premixed combustion 


15 






Fig. 8 Flame D (upper frame) and E (lower frame) thermo-chemical distributions taken at, xjdref — 1, 
xjdref = 15, and, xjdref — 30 from top to bottom. Model A, shown only in the first section corresponds 
to the red-dashed line; Model D, black-solid line; Symbols, experimental data. 


els C and D provide results in acceptable agreement with the experimental data, model D 
being slightly more accurate. Concerning the distribution of Z, an overall good agreement is 
achieved between the experimental data and the numerical results. The results obtained are 
in agreement with the results obtained using several combustion models using either RANS 
or LES approaches | 


4.2 Mixture fraction-conditioned results 

For an assessment of the combustion-model performance, a flow-field-independent compari¬ 
son with the experimental data 1341 would be needed. This can be obtained approximatively 
by analysing the Favre-averaged mixture fraction-conditioned thermo-chemical quantities 
obtained for both flames. Figure shows the conditioned temperature, methane, and water 
mass fractions distributions along three radial sections taken at dref, I5dref and 30dref, for 
flame D. Model D provides results in very good agreement with the experimental data El 
and always in better agreement with respect to the other three models. Focusing on the first 
section, one can see the remarkable differences between model A and model D; model A 
overestimates the methane consumption and, consequently, the temperature and the water 
productions. The differences fade through the three sections: the maximum temperature dif¬ 
ference is about 395^ at Z = 0.85 in the first section, whereas, it is about 137^ at Z = 0.215 
in the third section. The same conclusions are obtained analysing figure [^showing the Z- 
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Fig. 9 Flame D thermo-chemical Z-conditioned distributions taken at, xjdref = 1, xjdref = 15, and xjdref = 
30 respectively. Model A, red-dashed line; Model B, dashed-dotted line; Model C, dashed-dotted-dotted line; 
Model D, black-solid line; Symbols, experimental data. 



Fig. 10 Flame E thermo-chemical Z-conditioned distributions taken at, xjdref = xjdref = 15, mdxjdref = 
30 respectively. Model A, red-dashed line; Model B, dashed-dotted line; Model C, dashed-dotted-dotted line; 
Model D, black-solid line; Symbols, experimental data. 












































An SMLD joint PDF model for turbulent non-premixed combustion 


17 




Fig. 11 Flame D chemical species mass fractions Z-conditioned distributions taken at, xjdref — 1, x/dref = 
15, and xjdref = 30 respectively. Model A, red-dashed line; Model B, dashed-dotted line; Model C, dashed- 
dotted-dotted line; Model D, black-solid line; Symbols, experimental data. 


conditioned quantities at the same three sections for the flame E. It is noteworthy that all of 
the models overestimate the fuel consumption in the fuel rich side, even if the results are 
more accurate for model D. Such a behaviour is probably due to the model construction ne¬ 
glecting the interaction between turbulent mixing, chemistry, and radiative heat transfer CD. 
Figures[^and[^show the H 2 , CO and CO 2 Z-conditioned mass fractions, for flames D and 
E, respectively. Remember that the sum of the mass fractions of the three species together 
with the water mass fraction define the progress variable (see equation 0)- The figures are 
provided to check the model consistency with the physical phenomenon observed; in fact, 
one can see that in the farthest sections, where the reactions are fading (the progress variable 
tends to zero), the differences between the models are small with respect to the differences 
observed in the sections close to the burner. Einally, this analysis shows that all of th models 
provide quite good capability in the lean-fuel region of the flame, Z < 0.3. In the fuel-rich 
side, model D gives improved results, while slightly overestimate the temperature and the 
water productions due to the over prediction of the methane consumption, whereas model A 
provides unsatisfactory results. In this context it is remarked that even if model D is not 
able to replicate the accuracy of transported-PDE models, see e.g. I42l . its results appear to 
provide a remarkable improvement with respect to other EPV models. 


4.3 About the joint distribution of Z and A 

In this section we analyse the influence of two widely used hypotheses for the EPV models, 
namely, the statistical independence of Z and A and the j8-distribution PDE for P(Z). To this 
purpose, two points have been selected in the flow field of the Sandia flame E computed with 
model D and the corresponding values of the mean and the variance of both Z and C have 
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Fig. 12 Flame E chemical species mass fractions Z-conditioned distributions taken at, xjdref — 1, xjdref = 
15, and xjdref = 30 respectively. Model A, red-dashed line; Model B, dashed-dotted line; Model C, dashed- 
dotted-dotted line; Model D, black-solid line; Symbols, experimental data. 




Fig. 13 Distribution of Z and A at the point {xjdref ^yjdref) = (0; !)• solid line is the Psml,2 distribution 
and the dashed one the j8-distribution. 


been used to evaluate the probability distributions. For the first point, with normalized co¬ 
ordinates (xjdref .y!dref) — (0,1) (taken on the burner), the following values are obtained: 
Z = 0.2700, Z"3 = 0.0034, A = 0.9618 and A"^ = 0.0009. Note that the computational do¬ 
main, shown in figurej^ includes a part of the burner so that the flow properties at the section 
x/dref — 0 are computed whereas the boundary conditions are imposed at the inlet section at 
xjdref = —1.5. The resulting PDFs are shown in figure pA| It appears that the j8-distribution 
and the SMLD distribution of Z are in very good agreement. On the other hand, concerning 

























An SMLD joint PDF model for turbulent non-premixed combustion 


19 




(b) ix/dref,y/dref) = {20.^5,!) 


Fig. 14 Scatter-plot of the joint PDF with statistical independence hypothesis,P sml ,2 (2)P5ML ,2 (^), versus 
the joint PDF, Psml, 2 {'Z‘,A), at two different point in the Flame E flow fleld. The solid line is the bisector. 




Fig. 15 Distribution of Z and A at the point (xldref,y/dref) = (20.85, 1). The solid line is the Psml ,2 distri¬ 
bution and the dashed one the j8-distribution. 


the distribution of A, one can see that the two PDFs are quite different, providing different 
locations of the maxima and thus resulting in different evaluations of the mean values of 
the thermo-chemical quantities. It is worth noting that since j8(Z) and PsuL^i^) almost 
coincident, model C and model D differ only for the statistical independence hypothesis of 
Z and A. This issue is further analysed in figure(a), showing a scatter-plot corresponding 
to values for the mean and variance of Z and A computed at {x/dref,y/dref) — (0,1). The 
joint PDF is reported on both axes: the abscissa is evaluated using equation ( [23] ) whereas 
the ordinate is computed by the Bayes theorem with the statistical independence hypothesis. 
The bisector is also reported as reference, so that, if the statistical independence hypothe¬ 
sis were correct, one would obtain a cloud of points aligned with the bisector. Figure [T^a) 
clearly shows that for Flame E, the independence hypothesis is not appropriate in the region 
close to the burner. 
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Fig. 16 Distribution of Z and A at the point {x/dref^y/dref) = (1,1.21). The solid line is the Psml ,2 distribu¬ 
tion and the dashed one the j8-distribution. 


Consider, now, the second point (xjdref.yldref) — (20.85,1), far away from the burner. 
Ifee, the value of mean and variance are: Z = 0.3914, = 0.0397, A = 0.2074 and 

A"^ = 0.0822. At this location, the hypothesis that Z is distributed according to a j8-function 
results to be slightly less accurate, as shown in figure p~5] (left panel). The statistical inde¬ 
pendence hypothesis is slightly more appropriate than in the previous case; in fact one can 
see from figure [T^(b) that a non negligible part of the points are located near the bisector. 
Anyway, it still appears an incorrect hypothesis that should be abandoned in order to have 
an improvement in the combustion prediction. 

Now let us analyse the suitability of the assumption of the j8-distribution for Z and for A 
with an a-posteriori analysis. For this, three points are chosen diixjdyef — 1 (see the bottom- 
left panel of figure |^, where the difference between the results provided by model A and 
model D are found significant and the corresponding values of Z, Z"^, A, and are taken. 
The first point, y/dyef — 0.03, is correctly evaluated by both models and the values of mean 

and variance obtained using model D are: Z = 1, Z"^ = 0, A = 6.055 x 10“^, and A"^ = 0. It 
is clear that for these values j8(Z;Z,Z"2) collapses in 3{Z — Z) and j8(A;A,A"2) collapses 
in 5(A — A). In other words, the information related to the two values of the variance are 
lost, since the Dirac function is a one-parameter distribution depending only on the mean 
value. Using the and in the same way, due to the low 

values of the variances, the distributions collapse in two Dirac distributions. Very similar 
results are obtained for yjdref — 1-62, where the values of mean and variance found by 
model D are: Z = 0, Z^2 ^ 0, A = 0.001, and A^ = 6.449 x 10“^. 

The last point, yjdref = 1.21, is very interesting because of the different results obtained with 
the two models. Here the values of mean and variance found by model D are: Z = 0.221, 
Z"2 = 0.016, A = 0.557, and A"^ = 0.002. Figure [l^(le ft panel) shows that j8 (Z;Z,Z"^) and 

PsML,2{^\^,^"^) are very similar but PsAm2{M^,^"^) is more smooth than j8(A;A,A"^) 
(right panel). In this case the differences between the PDFs are such that ZuodeiD — 0.221 
and ZModeiA = 0.169 whit a relative error of about 23.5%. 
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5 Conclusions 

This paper analyses the state-of-the-art constitutive hypotheses usually adopted for the def¬ 
inition of the presumed probability density function (PDF) in flamelet progress variable 
(FPV) models, discussing their adequateness and feasibility. Four combustion models are 
considered, based on different choices for the PDF of the mixture fraction and the progress 
parameter. Among these models, a more general model is proposed evaluating the most 
probable joint distribution of the mixture fraction and the progress parameter using a closure 
technique based on an analytical evaluation of the Lagrange’s multiplier. The performance 
of the combustion models corresponding to different choices of the PDF are analysed by 
computing the Sandia flames D and E. The numerical results are also employed to study the 
validity of the hypothesis of statistical independence between the two variables. The results 
show that for RANS applications the proposed model provides a slight improvement in the 
prediction of the averaged thermo-chemical variables and a significant improvement in the 
evaluation of conditional quantities. In particular, the temperature distributions and its peak 
values are closer to the experimental data. This is also confirmed by an a-posteriori analysis 
showing that the introduction of joint PDF for Z and A is effective. The proposed model 
relies on a more robust theoretical basis, with a substantially unchanged computational cost. 
Finally, we can reasonably expect that the differences found with the proposed model in 
evaluating the conditioned quantity should be less pronounced when using LES. However, 
since the computational overhead is marginal it is recommended to consider the joint PDF 
also for LES applications. Future work will be devoted to study the influence of the proposed 
FPV model in the LES framework. 
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